package com.numericalmethod.algoquant.model.ralph2009.draft;

import java.util.Random;

/**
 * Trial implementation of 
 * 
 * "A Simulation Approach to Dynamic Portfolio
 * Choice with an Applicaion to Learning About Return Predictability"
 * (Brandt et al, 2005)
 * 
 * with one risky asset and no state variable.
 * 
 */
public class OneAssetSimulationTrial1_1 {

	// helper functions
	private static double average(double[] a) {
		double result = 0.0;
		for (int i=0;i<a.length;i++) 
			result+=a[i];
		result /= ((double) a.length);
		return result;
	}
	
	private static double stdev(double[] a) {
		double avg = average(a);
		double sum = 0.0;
		for (int i=0;i<a.length;i++) 
			sum+=Math.pow(a[i]-avg, 2.0);
		sum /= ((double) (a.length-1));
		return Math.sqrt(sum);
	}
	
	// CRRA Utility
	private static double utilityFunction(double w, double gamma) {
		
		return Math.pow(w, 1.0-gamma) / (1.0-gamma);
	}
	
	// CRRA differentiated once
	private static double utilityFunctionD1(double w, double gamma) {
		
		return Math.pow(w, -gamma);
	}
	
	// CRRA differentiated twice
	private static double utilityFunctionD2(double w, double gamma) {
		return -gamma * Math.pow(w, -gamma-1.0);
	}
	
	// (TESTING DISPLAY)
	private static void printComparison(double[] nextExcessReturn, double estimatedOptimalWeight, 
			double mu, double sigma, double riskFree, double dt, double gamma) {

		System.out.println("==== COMPARISON ====");

		double averageExcessReturn = average(nextExcessReturn);
		double returnVolatility = stdev(nextExcessReturn);
		System.out.println("Next Step Excess Return = "+averageExcessReturn);
		System.out.println("Next Step Vol = "+returnVolatility);
		System.out.println();
		
		// close form for Mean/Variance portfolio optimization on CRRA Utility 
		// ???CORRECT????
		double closeFormWeight = 1.0/gamma * ((mu-riskFree)/(sigma*sigma));
		
		// (TESTING DISPLAY)
		// the portfolio return mean and vol using our estimated weights
		double ptfMean = averageExcessReturn*estimatedOptimalWeight+(1.0-estimatedOptimalWeight)*riskFree*dt;
		double ptfVol = returnVolatility * estimatedOptimalWeight;
		System.out.println("Using Simulation Estimated Weight: "+estimatedOptimalWeight);
		System.out.println("portfolio mean = "+ptfMean);
		System.out.println("portfolio vol = "+ptfVol);
		System.out.println("portfolio mean/variance = "+(ptfMean/(ptfVol*ptfVol)));
		System.out.println();

		// the portfolio return mean and vol using the close form weights.
		ptfMean = averageExcessReturn * closeFormWeight+(1.0-closeFormWeight)*riskFree*dt;
		ptfVol = returnVolatility * closeFormWeight;
		System.out.println("Using Mean/Variance Analytic Weight: "+closeFormWeight);
		System.out.println("portfolio mean = "+ptfMean);
		System.out.println("portfolio vol = "+ptfVol);
		System.out.println("portfolio mean/variance = "+(ptfMean/(ptfVol*ptfVol)));
		System.out.println("====================");
		System.out.println();
		System.out.println();
		
	}
	
	private double[] optimalWeight;
	private double analyticWeight;
	
	public OneAssetSimulationTrial1_1() {
	}
	
	public double[] getOptimalWeight() {
		return optimalWeight;
	}
	
	public double getAnalyticWeight() {
		return analyticWeight;
	}
	
	
	/**
	 * Solve for optiomal weights
	 * 
	 * @param rand random generator to be used
	 * @param mu annualized return rate
	 * @param sigma annualized volatility
	 * @param riskFree risk free rate
	 * @param pathCounts simulation counts
	 * @param stepsCount the algorithm requires at least 2 steps
	 * @param gamma CRRA utility parameter
	 * @param dt time step size
	 */
	public void solve(
			Random rand,
			double mu,
			double sigma,
			double riskFree,
			int pathCounts,
			int stepsCount,
			double gamma,
			double dt) {			
		
		// close form solution
		analyticWeight = 1.0/gamma * ((mu-riskFree)/(sigma*sigma));
	
		// unknowns to be solved by simulations
		optimalWeight = new double [stepsCount-1];
		
		// generate the excess returns
		double[][] excessReturns=new double [pathCounts][stepsCount];
		for (int m=0;m<pathCounts;m++) {
			for (int t=0;t<stepsCount;t++) { 
				excessReturns[m][t] = (mu )*dt+sigma*Math.sqrt(dt)*rand.nextGaussian() - riskFree * dt;
			}
		}
		
		// move backward in time
		for (int t=stepsCount-2;t>=0;t--) {
			
			double[] aTPlus1 = new double [pathCounts];
			double[] bTPlus1 = new double [pathCounts];
			
			// ??? w(t) = 1.0 for all??? as mentioned at the bottom of page 8 in the paper
			// because we are using CRRA?
			double currentWealth = 1.0;
			// If not, I don't know how to compute the wealth at time t 
			// because it requires the current time "optimal weight" which is not yet computed
			// (budget constraint mentioned in equation (2) in the paper)

			for (int m=0;m<pathCounts;m++) {
				// some factor variables that reuse again and again
				double factor1 = 1.0;
				double factor2 = 1.0;
				for (int s=t+1;s<stepsCount-1;s++) {
					double periodReturn = (1.0 + optimalWeight[s] * excessReturns[m][s+1] + riskFree * dt); 
					factor1 *= periodReturn;
					factor2 *= Math.pow(periodReturn, 2.0);
				}
				
				// (equation 15)
				double estimatedTerminalWealth = currentWealth * (1.0 +riskFree * dt) ;//* factor1;
				
				
				// (equation 18)
				aTPlus1[m] = utilityFunctionD1(estimatedTerminalWealth, gamma) * factor1 * (excessReturns[m][t+1]);
				bTPlus1[m] = utilityFunctionD2(estimatedTerminalWealth, gamma) * factor2 * Math.pow(excessReturns[m][t+1], 2.0);

				
			}
			
			
			// (equation 19 and 20)
			// we don't have any state variable for simplicity
			// the cross paths regression becomes just plain average of all generated values at the same timestep
			double expectedATPlus1 = average(aTPlus1);
			double expectedBTPlus1 = average(bTPlus1);
			

			// the optimal weight (equation 21)
			optimalWeight[t] = (-1.0 / ( expectedBTPlus1 * currentWealth)) * expectedATPlus1;
			
			// (TESTING DISPLAY)
			System.out.println("At time step "+t+": Estimated Optimal Weight = "+optimalWeight[t]);
			double[] nextExcessReturn = new double [pathCounts];
			for (int i=0;i<pathCounts;i++) {
				nextExcessReturn[i]=excessReturns[i][t+1];
			}
			printComparison(nextExcessReturn, optimalWeight[t], mu, sigma, riskFree, dt, gamma);
			// (TESTING DISPLAY END)

		}
	}
	

	public static void main(String[] args) throws Exception {
		
		// fixed seed for comparison purpose
		Random rand=new Random(2965179875022589598L);

		// risky asset SDE parameter.
		// r = mu dt + sigma dW
		// annualized return rate
		double mu = 0.92/100;
		// annualized return volatility
		double sigma = 5.24/100;
		
		// annualized risk free rate
		double riskFree =0.005;

		// simulation counts
		int pathCounts=100000;
		// the algorithm requires at least 2 steps
		int stepsCount=10;
		
		// CRRA utility parameter
		double gamma = 5.0;

		// time step is 1 business day
		double dt=1.0/250.0;


		OneAssetSimulationTrial1_1 t=new OneAssetSimulationTrial1_1();
		t.solve(rand, mu, sigma, riskFree, pathCounts, stepsCount, gamma, dt);
		
		
		System.out.println("Simulation Finished");
	}
	
	
	
	
}
